home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / MFLOAT10.ZIP / BESSEL.PAS < prev    next >
Pascal/Delphi Source File  |  1993-04-28  |  2KB  |  84 lines
  1. PROGRAM bessel(input, output);
  2. { *** This test program shows the advantages of mfloat numbers. *** }
  3. { *** The bessel function J0(x) is calculated by the series.    *** }
  4.  
  5. USES pfloat, crt;
  6.  
  7. {-------------------------------------------------------------------------}
  8.  
  9. FUNCTION J0(x : extended) : extended;
  10.  
  11. CONST epsi = 1e-20;
  12.  
  13. VAR sum, sqrx, prod, mepsi : mfloat;
  14.     i                      : integer;
  15.  
  16. BEGIN
  17.   ldtomf(mepsi, epsi);
  18.   ldtomf(sqrx, x);
  19.   ldexpm(sqrx, -1);
  20.   sqrm(sqrx);
  21.   Getonem(sum);
  22.   Getonem(prod);
  23.   i := 0;
  24.   REPEAT
  25.     i := i + 1;
  26.     multm(prod, sqrx);
  27.     divi(prod, i);
  28.     divi(prod, i);
  29.     IF odd(i) THEN
  30.       subm(sum, prod)
  31.     ELSE
  32.       addm(sum, prod);
  33.   UNTIL gtm(mepsi, prod);
  34.   J0 := mftold(sum);
  35. END;
  36.  
  37. {-------------------------------------------------------------------------}
  38.  
  39. FUNCTION J0x(x : extended) : extended;
  40.  
  41. CONST epsi = 1e-20;
  42.  
  43. VAR sum, sqrx, prod : extended;
  44.     i               : integer;
  45.  
  46. BEGIN
  47.   sqrx := sqr(x / 2);
  48.   sum := 1;
  49.   prod := 1;
  50.   i := 0;
  51.   REPEAT
  52.     i := i + 1;
  53.     prod := prod * sqrx / i / i;
  54.     IF odd(i) THEN
  55.       sum := sum - prod
  56.     ELSE
  57.       sum := sum + prod
  58.   UNTIL epsi > prod;
  59.   J0x := sum;
  60. END;
  61.  
  62. {-------------------------------------------------------------------------}
  63.  
  64. VAR accuracy : integer;
  65.     x        : extended;
  66.  
  67. BEGIN
  68.   ClrScr;
  69.   writeln('              Calculation of the bessel function J0(x)');
  70.   writeln('              ========================================');
  71.   writeln;
  72.   x := 100;
  73.   writeln('mantissa words        result     (x = ', x:5:1, ')');
  74.   writeln;
  75.   FOR accuracy := 1 TO 15 DO
  76.   BEGIN
  77.     Setmantissawords(accuracy);
  78.     writeln(accuracy:3, '                J0(',x:5:1,') = ', J0(x));
  79.   END;
  80.   writeln;
  81.   writeln('IEEE arithmetic    J0(',x:5:1,') = ', J0x(x));
  82.   writeln;
  83.   writeln('The accuracy of 12 mantissa words is needed to get a good result!');
  84. END.